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THEORETICAL STUDIES OF SOLAR LASERS AND CONVERTERS 


By 

John H. Heinbockel* 

Research for the period consisted of developing and refining the 
continuous flow laser model program. In addition, a working model for a two 
pass continuous wave amplifier was developed. The following is a summary of 
the mathematical development of a two pass amplifier for the n-C 3 FyI iodine 
laser. The geometry of the amplifier is illustrated in Figure 1. In this 

- power density into the amplifier (watts/cm^), 

- transmission coefficient, 

• radiation energy from I* (watt-sec), 

- speed of light, 

- output power density, 

- reflection coefficient for mirror, 

- photon density for wave motion in positive direction (cm J ) , 

- photon density for wave motion in negative direction (cm"- 5 ). 

The input power density P^ n is assumed to be known. The photon density 
satisfies for all values Z 

P+(Z) p- (Z) — Kq 
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where Kq is a constant. At Z - 0 we have 


P+(0)p-(0) - K 0 
and at Z - L we have 

P+(L) p-(L) - Kq. < 3) 

The value of p+(0) is determined from the assumed value for the input 
power density P^ n and from this we calculate 


pH 0) - 


P. W 
in 

£ C 


and consequently from the relation (2) we have 


(4) 


pH 0 ) 

In the above relations W is the transmission coefficient of the Brewster 
windows and K 0 is an assumed initial value . At Z - L we have R 2 - W 2 R L as 
the reflection coefficient from the Brewster window. We then are able to 
calculate the return photon density from the Brewster window at Z - L. We 
find that this density should have the value 
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■V 


p-( L) - R 2 p+(L) . 

Combining the equations (2) (3) and (5) we determine that 


R 2 P+ 2 (L) - K 0 


or 


pH L) 



P+(0)p-(0) 


( 6 ) 


Using the value (6) as the theoretical value for p+(L) at Z - L we can com- 
pare this value with the calculate value for p+(L) obtained from an integra- 
tion of the differential equations defining the chemical reactions occurring 
in the amplifier. We adjust the value for p-(0) until 




w 


theoretical calculated 


< E 


(7) 


where E is an error criteria. 

The above procedure is a shooting method for the determination of p-(0). 
When the error criteria is satisfied, the output power is obtainable from the 


relation 


P - p-( 0) W e C 
out 

and the gain is calculated from the relation 


gain - 10 log 1Q 



(decibels) . 


(9) 
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The chemical reactions occurring in the laser tube are listed along with 
the values used for the rate constant. The assumed reactions are: 
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* 
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12 
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* 

I+I 

(51% 

of 

the 

time) 

12 

+ 

hw 
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(49% 

of 

the 
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Two Body Reactions 

K, 


* 

I +R 


I+R 

R+R 

R+RI 

r+i 2 

R+RI 

* 

I +RI 


-> RI 




K. 


K, 


RI 


r 2 +i 

RI+I 


r 2 +i 

R+I„ 


I+RI 


* X 1 
I +RI — 


★ 

I +1, 


* Q 3 

I +R 

* ^4 

I +R„ 


* 

I +1 


i 2 +r 

I+RI 

1+1 2 
I+R 

i+r 2 

I+I 



-14 

- 1.0 (10 V ) 

- 0.8 ( 10 ' 11 ) 

-12 

- 2.0 (10 ) 

- 3.0 (10' 16 ) 

- 1.0 (10' 11 ) 

- 0.0 

- 3.0(10' 9 ) 

-23 

- 1.6(10 

- 1.7(10' 17 ) 

- 3 .8(10- 11 )e- 4 - 4(10 ' 3><T - 300> 

- 3.7(10' 18 ) 

- 4.7(10' 16 ) 

- 1 . 6 ( 10 ' 14 ) 


Three Body Reactions 

I*+I+RI I„+RI 

C 2 

I+I+RI — — » I 2 +RI 


. C x - 1.6(10' 33 ) 

, C 2 - 5.7(0- 3 V 31 °/ r 
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* 

I +I+I, 


I+I+I, 


I+I+R 

* C 6 
I+I +r 2 


I 2 +I 2 

, c 3 - 0.0 

V X 2 

■ c 4 - io a 

I 2 +R 2 

, C 5 - 8.0(10 

I 2 +R 2 

, c 6 - 0.0 


-33, 1310/r 


where 


a - -29.439 -5.844 log 


10 


300 


-2.163 


log 


10 


300 


Let B denote a scale factor and let y^, i— 1,7 denote the normalized 
values of concentration. Define the quantities: 

Xi - Byi - [RI] 

X 2 - By 2 “ [R] 

X 3 - By 3 - [R2 1 

X 4 - By 4 - [I 2 ] 

X 5 - By 5 - [I*] 

X6 - By 6 - [*] 

X 7 - By 7 - p+( Z) 

X 8 - p - (Z)— Ko/Xy 
p - p(Z) - X 7 +X 8 


where [ ] denotes concentration of substance, R is the iodide perflouride 
radical n-C 3 F 7 , I is iodine and I* denotes the excited state of iodine. For 
w the flow rate of n-C 3 F 7 I in the Z direction (cm/sec) and p the total photon 
density, the differential equations defining the chemical kinetics of the 
amplifier are given by: 
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dy^ 

dz 

^2 

dz 


K l B y2 y 5 +K 2 By 2 y 6^1 y r K 4 B V2 +K 5 By 2 y 4- K 7 By 5 y r K 6 By 2 y r K 8 By l y 6- y l- 
€ l y l * K l By 2 y 5 * K 2 By 2 y 6 ' 2K 3 By 2 ’ K 4 By l y 2 ' K 6 By l y 2 - K 5 By 2 y 4 


+ K 7 By 1 y 5+ K 8 By 1 y 6 -y 2 - y 2 /r 2 

dz 


dy 3 2 dw 

« - K 3 By 2 + K 6 B yi y 2 + ^By^^ - - y 3 A 3 

dz dz 


dz 


- C l B2 5'iy5J , 6 tC 2 B V6 2+C 3 B V55 r 6- K: 4 B2j, 4 y 6 
• S By 2 y 4 +C 5 B2y 6 2y 3 +K 8 By 6 y l +C 6 B2y 6 y 5 y 3- y 4 


2 -e 2 V K 7 By i y 5 

dw . 

— - V'4 


• — ' « 1 y 1 *°-5H 2 y 4 -K l By 2 y 5 -C 1 B 2 y 1 y 5 y 6 -C 3 B 2 y 4 y 5 y 6 -Q 1 By 1 y 5 -Q 2 By 4 y 5 

dz 

-Ca.p.(y 5 - 1 y 6 ) + K 6 By 2 y 1 -Q 3 By 5 y 2 -Q 4 By 5 y 3 -Q 5 By 5 y 6 

2 du , 

-K 7 By 5yi -c 6 B y 6 y 5 y 3 -y 5 — - y 5 A 5 

dz 

» — - l-« « 2 V ( ’l y l y 5 +<5 2 By 4 y 5- 2C 5 B2y 6 2y 3- K 8 By 6 y l + ^ <y 5‘ ; V 
dz 1 

• C l B2y l y 5 y 6' 2C 2 B2y l y 6 2 ' C 3 B2y 4 y 5 y 6' 2C 4 B2y 4 y 6 2 ' K 2 By 2 y 6 +K 4 By l y 2 
+<) 3 By 5 y 2 4<5 4 By 5 y 3 +<, 5 By 5 y 6 +K: 5 By 2 y 4- C 6 B2y 6 y 5 y 3- y 6 ^ ’ V r 6 
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dy 7 1 

— - Vl (y 5' - y 6 )B< ’ 

dz l 

Normalizing the length of the amplifier by using the substitution Z - 
SL, 0 < S < 1 the differential equations are solved numerically using a 
Runge-Kutta integration scheme. Appendix A contains the computer code for 
the simulation of a two pass continuous wave amplifier. Appendix B contains 
a representation of the output from the code. This Appendix contains figures 
illustrating how concentrations change with distance. Current research is 
directed toward the refining of this program and performing a systems study 
with the parameters involved in the code construction. 
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APPENDIX A 


COMPUTER PROGRAM FOR SIMULATION OF TWO PASS CONTINUOUS WAVE AMPLIFIER 


PROGRAM TPAMP2 ( INPUT, OUTPUT, TAPE5= INPUT, TAPE6=0UTPUT, TAPE8) 

C 

C 

C MAIN PROGRAM 

C TWO PASS CONTINUOUS WAVE AMPLIFIER 

C 

C AMPPAR IS THE DATA FILE WITH NAMELIST VARIABLES DEFINED 

C 

C COMPRESSIBLE FLOW LASER MODEL WHICH INCLUDES 

C EQUATION OF STATE 

C CONTINUITY EQUATION 

C MOMENTUM EQUATION 

C ENERGY EQUATION 

C 

C TAPES IS FOR OUPUT DATA RENAME AND SAVE TAPES IF YOU WANT 

C TO SAVE THE OUTPUT DATA 

C PROBLEM IS SCALED IN THE AXIAL. DIRECTION LENGTH IS LC 

C SCALED LENGTH IS ZERO TO ONE 

C 

C SPECIFIC HEAT AT CONSTANT VOLUME IS FUNCTION OF TEMPERATURE 

C OBTAINED FROM LEAST SQUARES FIT OF EMPIRICAL DATA 

C FOR N-C3F7I IODINE LASER 

C 

C USED IN PREDICTION MODE GIVEN CONDITIONS AT Z=0. 

C 

c 

C THIS VERSION CONTAINS AN AUTOMATICALLY CHOSEN VARIABLE STEP 

C SIZE. IT ALSO CONTAINS WALL. EFFECT REACTIONS. (NOT USED THOUGH) 

C 

C EFFECTS OF TEMPERATURE, FLOWRATE, DENSITY AND PRESSURE VARIATIONS 

C ARE CONSIDERED IN SUBROUTINES ARREN, PFLOW AND FLOW 

C 
C 

COMMON / BLK2 / K 1 , K2 , K3 , K4 , K5 , K6 ,K7, K8 ,C1,C2 , C3 , C4 , C5 
COMMON /'BLK27 /Q1 , Q2 , Q3, Q4 , Q5 

COMMON /BLK3/B, B2, B3, C, AOO, BOO, EPSNU, OMEGA, C6 
COMMON/ BLK4/CHS 110, CHS 1 20 , ABARO , Z 1 BAR , LC , XKO 
COMMON / BLK7 / ABC , COO , CO , OMEG 1 , P , R 1 , R2 , TM , XNRHO 
COMMON/BLK8/ZOL, ZE, NG, TO, RAD, A, PIN, W 

COMMON /BLK 1 0/CF 1 , CF2, CF4, QFO, RSTAR, ZL, L, SF1 , SF2, AR, AAO, BBO 
COMMON /BLK22/ AD, VI , V2, GG 
COMMON /BLK23/ WO, ETAO, PTO, FRAC 
C0MM0N/BLK28/TAU2, TAU4, TAU5, TAU&, TAU3, SIG 
COMMON /BLK299/TTT2, TTT3, TTT4 , TTT5 , TTT6 
REAL K 1 , K2 , K3 , K4 , K5 , K6 , K7 , K8 
REAL LC, L 
C 
C 

WRITE (6, 123) 

123 FORMAT ( IX, 20H START OF PROGRAM ) 
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DEFAULT VALUES 
IEND=0 
I COUNT =0 
NG=0 

FRAC=U. U 
AR=2. 

TLE-39 . 0 
PT 0=20 . 0 
W0= 15.8 
ETA0=. 33 
ZOL-7 . 50 
T 0—300 . 

AA0=1 47 , 23 
A AO IN J/K MOL 
BBC IN 1/K 
BBC- . 00 1 2 

AAO, BBO USED IN LEAST SQ FIT OF SP. HEAT CONSTANT VOLUME 
LC= 15.0 
L.C IN CM 
XNRHO= 1 , 0 
RAD=0. 0 
VI =0.0 
V2=0 . 0 
TTT2=1 . 0E18 
TTT3= 1 . OE 1 8 
TTT4= 1 . OE 1 8 
TTT5= 1 . OE 1 8 
TTT6=1 . OE 1 8 
C00=. 5E12 
R1=0. 0 
R2=. 975 
T 0=300 . 

P I N= 1 . 0 
W=0 . 98 
ZE=5 . 0 
A=2. 0 
TM= 1-R2 
OMEG 1=5000. 

0MEG1 IN CM/S 
C0N=2. 0E4 


ORIGINAL PAGE IS 
OF POOR QUALITY 


“ NAMEL I ST/PARAM / PTO , OMEG 1 , CON , COO , R 1 , R2 , TM , XNRHO , LC , ZOL , A , ZE , 

1 AR, I END, FRAC, TO, RAD, VI , V2, TTT2, TTT4, TTT5, TTT6, TTT3, W, PIN 
C 

c 

_ C ASSUME PRESSURE PTO (TORR) AND TEMPERATURE TO (DEG K) 

C ARE GIVEN FOR LEFT END. CALCULATE ETAO TO SATISFY GAS LAW 

55 CONTINUE 

IF ( IEND .EG!. 1) GO TO 600 

IF(NG.EQ.S) GO TO 600 
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ORIGINAL PAGE IS 
OF POOR QUALITY 


W--TRANSM I S3 1 ON COEFF I C I ENT 
P IN- INPUT POWER DENSITY (WATTS/CM**2> 

P -PRESSURE, TORR 
AF -LASER BEAM DIAMETER CM 
A-RADIUS OF TUBE (CM) 

RAD=RAD I US OF TUBE WHERE CALCULATIONS ARE DONE. 

T=TEMF'ER ATURE DEG K 
TO-TEMFERATURE LEFT END DEG K 
P“0= PRESSURE TORR LEFT END 
WO-FLOW RATE LEFT END M/SEC 
E T AO=DENS I TY OF GAS KG/M**3 

OMEG INFLOW RATE, CM/SEC , MAXIMUM FLOW RATE AT RAD=0 
TAU2j TAU4, TAU5, TAU6 ARE LIFETIMES FOR RADIAL DIFFUSION OF 
T T2, TTT3, TTT4, TTT5, TTT6 ARE LIFETIME VALUES FOR TAU 
THE QUANTITIES CR3 , C 12] , C 1*3 AND III. 

TAU3 IS DIFFUSION COEFFICIENT FOR CR2H 
CON=PEAK CONCENTRATION , SOLAR CONSTANTS 
COO- INITIAL GUESS AT RHO-F'LUS AT ZERO 
WHICH IS SQUARE OF (COO*Rl) 

Rl= REFLECTIVITY AT LEFT END 
R2= REFLECTIVITY AT RIGHT END 
ZE=DISTANCE FOR LIGHT INTENSITY TO DIMINISH BY FACTOR 1/E 
TM= TRANSMISSION COEFFICIENT (OUTPUT MIRROR) 

ZOL=POINT ALONG AXIS WHERE MAXIMUM ILLUMINATION OCCURS 
IN THE CASE ILLUMINATION IS A BELL SHAPED CURVE 
IN THE CASE OF A SQUARE WAVE, 2*Z0L IS CUT OFF POINT 

THE POINT 2*ZOL IS WHERE ILLUMINATION BEGINS TO DIMINISH 
LC=LENGTH OF LASER CAVITY IN CM 
S= SCALED LENGTH WHICH GOES FROM ZERO TO ONE 
.LENGTH IS SCALED FROM ZERO TO ONE 
ZL=2*Z0L LENGTH WHICH IS ILLUMINATED 
FRAC=FRACTION OF PEAK CONCENTRATION WHICH GOES 
INTO HEAT 

DEFAULT VALUE IS FRAC-0.024 I.E. 2.40 PERCENT 
H IS STEP SIZE FOR NUMERICAL INTEGRATION OVER SCALED LENGTH 


:>1 CONTINUE 
F-PTO 
T=TO 

ETAO=F'* ( 1 . 01325E5) *296. / (8314. 3*760. *T) 
FT 0 IS PRESSURE IN TORR 
ET AO IS DENSITY IN KG/M**3 
ETAO IN KG/M**3 


IC0UNT=IC0UNT+1 
READ (5, PAR AM) 

IF (EOF (5) ) 600, 601 
00 WRITE (6, 603) 

03 FORMAT ( IX, 28HEND OF FILE ENCOUNTERED-STOP ) 

CALL PSEUDO 
DO 10 JJJ=1,9 
CALL GRAPHS (JJJ) 

0 CONTINUE 

STOP 1313 
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c, 

c 


c 

c 

1 Ci 8 
19? 


C 

c 


c 

c 


c 

c 

c 


cc 

cc 

702 

C 


C 


701 


708 


ORIGINAL PAGE IS 
OF POOR QUALITY 


F I N TORR 
~ IN DEG K 
CMIN=1.0E4 
CMAX=1 . 0E20 
ZL=2*ZQL 
L=LC 

c :i=con 
c: i=con 

C-LL COEFFS 
.-.0 IN M/SEC 
;*- = 0 M E G A / 1 OO . 

Z- 4 EG1 AND OMEGA ARE IN CM /SEC 
^- r :ITE (6, 198:» 

F IF:MA T ( / / / ) 

W- ITE (6, 199) ZE, ZOl , CON, OMEGA, FR AC, C00,R1 , R2 
FORMAT (IX, T5,5HZE * , F10. 3, T30, SHZQL = , E15 

1 E 1 5 . 7 , T 90 , 8HQMEGA = ,E15„7,T1 04 , 7HFRAC = ’ , F9 

2 1 : , T5« tsHCOO - 


F,T 

7, T55 . 6HCC.N 
A .« / , 

, F 


E15.7, T30, 6H R1 - ,Fl0.7, T55,tH R2 

F L 0 . 3 ! 

EQUATIONS 


31, 4N F’ — , E 1 5 • 7 , T 1 05 , /H I EMF = , F L 0 . 3 


SE“ UP COEFFICIENTS IN DIFFERENTIA 
SET PRINTER OFF 
I PR I NT =0 PRINTER OFF 

I PRINT = 1 PRINTER ON 

I =, RINT=0 

SET STEP SIZE H=.025 UNITS 
►- =*0 • 025 

INTEGRATE DIFFERENTIAL EQUATIONS FROM 3=0 TO S=1.0 
X 1 =000 

CALL I NTEG ( I PR I NT , H ) 


INTERVAL HALVING SCHEME 


Y 1 =ABC 

IF(Yl.LT.O) PER-5.0 
IF(Yl.LT.O) PER- 1.1 
IF(Yl.GT.O) PER*. 1 
IF (Y 1 „ GT . 0) F'ER= * 9 
CONTINUE 

IN THIS VERSION COO IS ASSUMED VALUE FOR RHO-(O) 

COO= (PER) *C00 

IF (COO . LT. CMIN) STOP 5432 
IF (COO .GT. CMAX ) STOP 2345 
X2=C00 

CALL I NTEG ( I PR I NT, H) 

Y2=ABC 

ABC IS XKO INITIAL MINUS XKO CALCULATED (PERCENT ERROR) 

IF( (Y1*Y2> .LT. 0) GO TO 701 

X1=C00 

Y 1 =Y2 

GO TO 702 

CONTINUE 

W 1 = . 6 

W2= . 4 

COO=W2*Xl+Wl*X2 

CONTINUE 

CALL I NTEG ( IPRINT, H) 
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704 

C 


705 

C 


ccrc 
•-J wJ 


C 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


X3=C00 

V3=ABC 

IF vABS(Y3> . LT.0.05) GO TO 555 
CONTINUE 

IF ( ( Y1*Y3) . LT . 0) GO TO 705 
Y1 Y3 ARE OF THE SAME SIGN 
X i=X3 

Y I =Y3 

C X)=.5*<X1+X2) 

GO TO 708 
CONTINUE 

Vi Y3 ARE OF OPPOSITE SIGN 
X D=X3 

Y 0— Y3 

C 30= . 5* ( X 1 + X 2 ) 

GO TO 70S 
I PR I NT - 1 
-=0. 025 

CALL I NTEG ( I PR I NT , H ) 

GO TO 55 
END 
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SUBROUTINE GRAPHS* JJJ) 


AS THE NAME SUGGESTS THIS SUBROUTINE PLOTS GRAPHS OF THE 
NUMERICAL OUTPUT 

MULTIPLE PLOTS ON EACH GRAPH UP TO MAXIMUM OF EIGHT 


COMMON/BLK4/CHSI 10, CHS I 20 , AB ARO , Z 1 BAR , LC , XKO 
COMMON /BLK8/Z0L , ZE, NG, TO, RAD, A, PIN, W 

COMM0N/BLK30/DATA (50, 50) , NDMAX , FLRATE (3) , FDPTDAT (50, 40) 
DIMENSION X ( 50 ) , Y C 50 ) , Y Y ( 50 , 8 ) 

REAL LC 


JJJ=1 , 9 
JJJ=1 PLOT 
JJJ=2 PLOT 
JJJ=3 PLOT 
J J J=4 PLOT 
JJJ=5 PLOT 
J J J=6 PLOT 
J J J=7 PLOT 
JJJ =8 PLOT 
J J J=9 PLOT 


R VS Z 
12 VS Z 
I * VS z 
I vs z 

INV VS Z 
FLOWRATE VS Z 
DENSITY VS Z 
PRESSURE VS Z 
TEMPERATURE VS Z 


DATA ARRAY IS BY COLUMNS 

Z,R, 12, I*, I, INV, Z,R, 12, I*, I, INV, . . . 

ND=NUMBER OF DATA POINTS, =NUMBER OF ROWS IN DATA ARRAY 
NG=NUMBER OF CURVES PER GRAPH 
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C IC IS CODE TO DETERMINE NUMBER OR GRAPHS TO PLOT 

C IC=0 FOR MORE THAN ONE GRAPH 

C IC=1 FOR LAST GRAPH (USED FOR ONLY ONE GRAPH) 

C 

C 

IC=0 

NLAST =NG- 1 
SXC=5. 0 
SYC=5, 0 
3 VC I =7 . O 

C SXC=LENGTH OF XSCALE 

C SYC=LENGTH OF YSCALE 

C S V CI =LENGTH OF YSCALE FOR INVERSION 

IF (NLAST . EQ. 0) IC=1 
C PLOT JJJ VS Z 

IF ( J J J . GE. 6) GO TO 700 
DO 10 1 = 1 , NDMAX 
X ( I > =DAT A ( I , 1 ) 

DO 20 J=1,NG 
NN= < J-l ) *6 
NCOL=NN+ 1 + J J J 
YY ( I , J ) =DAT A(I, NCOL ) 

20 CONTINUE 

CO CONTINUE 

C 

C FIND YMAX , YMIN 

C 

C 

DMAX=0 . 0 
YMAX=0. 0 
YM I N=0 . 0 
ZMIN=0. 0 
C 

C FIND YMAX AND YMIN 

C 

DO 70 1=1, NDMAX 
DO 71 J= 1 , NG 

IF ( YY ( I , J ) .GT. YMAX) YMAX=YY(I,J) 

IF ( YY ( I , J ) .LT. YMIN) YMIN=YY(I,J) 

71 CONTINUE 

70 CONTINUE 

C 

C PLOT FIRST DATA CURVE 

C 

C JJJ=1 

ZMAX= 1 . 00 
DO 40 1=1, NDMAX 
Yd) =YY (1,1) 

40 CGNTINUE 

L0GMAX=1+AL0G10 (ABS (YMAX) ) 

POWER1=10. **LOGMAX 
TEST = YMAX /POWER 1 
YMAX X=0. 025 

IF (TEST .GT, 0.025) YMAXX=.05 
IF (TEST .GT. 0.05) YMAXX=0. 10 
IF (TEST .GT. 0.10) YMAXX=0.25 
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506 


507 

C 

C 


50 


C 

C 

C 


1 00 


C 

C 


L. 200 


c 

c 

c 


IF (TEST .GT. 0.25) YMAXX=0.50 
IF (TEST .GT. 0.5) YMAX X= 1.00 
YMAX»YMAXX*P0WER1 
IF ( J J J . NE . 5) GO TO 507 
CONTINUE 

IF (YMIN .GE. 0) GO TO 507 

IF (YMIN . LT. 0.0) L0GMAX2=1+AL0G10 (ABS ( YMIN) ) 

P0WER=10. **L0GMAX2 
TEST =ABS (YMIN) /POWER 
VMINN=0. 001 

I F ( TEST . GT . 0.00 1 ) YM I NN-0 . 005 
I F ( TEST . GT . 0 . 005 ) YM I NN=0 . 0 1 
IF (TEST .GT. 0.01) YMINN=0.05 
IF(TEST .GT. 0.05) YMINN=0.1 
IF (TEST . GT. 0.10) YMINN=0.5 
I F ( TEST . GT . 0.5) YM I NN= 1 . 0 
YM I N=-YM I NN*POWER 
SYC=SYCI 
CONTINUE 
Z M I N=0 . 0 

IF(JJJ.GT.l) GO TO 50 
JJJ=1 PLOT R VS Z 

CALL INFOPLT ( IC, NDMAX , X , 1 , Y, 1 , ZMIN, ZMAX , YMIN, YMAX , 1.0, 

1 22, 22HZ , SCALE DISTANCE, 

2 11, 1 1 H CC3F73 ,0, 

3 SXC, SYC, 0. 75, 0. 75) 

GO TO 600 
CONTINUE 

IF ( J J J . GT. 2) GO TO 100 
J J J=2 PLOT 12 VS Z 

CALL INFOPLT ( IC, NDMAX, X, 1,Y, 1 , ZMIN, ZMAX , YMIN. YMAX , 1.0, 

1 23, 23HZ , SCALE DISTANCE, 

2 9 , 9H C 123 ,0, 

3 SXC, SYC, 0.75, 0.75) 

GO TO 600 
CONTINUE 

IF ( JJJ . GT. 3) GO TO 200 
J JJ=3 PLOT I* VS Z 

CALL INFOPLT (IC, NDMAX, X, i,Y, 1 , ZMIN, ZMAX, YMIN, YMAX, 1.0, 

1 23 , 23HZ , SCALE DISTANCE, 

2 9, 9H C 1*3 ,0, 

3 SXC, SYC, 0.75, 0.75) 

GO TO 600 
CONTINUE 

IF (JJJ . GT. 4) GO TO 300 
J J J=4 PLOT I VS Z 

CALL INFOPLT ( IC, NDMAX , X , 1, Y, 1, ZMIN, ZMAX , YMIN, YMAX , 1.0, 

1 23, 23HZ , SCALE DISTANCE, 

2 8, 8H C 1 3 ,0, 


v* 
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c 

c 

c 


600 

C 

c 


501 


500 
500 1 
C 

c 

c 


60 


601 


700 

C 


702 

C 

C 

701 

C 


704 

703 


3 SXC,SYC, 0.75, 0.75) 
GO TO 600 
CONTINUE 


J J J=5 PLOT INV VS Z 

CALL INFOPLT ( IC, NDMAX, X , 1, Y, 1, ZMIN, ZMAX , YMIN, YMAX , 1.0, 

1 23, 23HZ , SCALE DISTANCE, , 

2 24 , 24H Cl*]-. 5*C 1 3 ,0, 

3 SXC,SYC, 0.75, 0.75) 

CONTINUE 

PLOT REST OF CURVES OR EXIT IF ONLY ONE CURVE 
NLAST =NG- 1 

I F (NG ,EQ. 2) GO TO 5001 
IF (NLAST .EQ. 0) GO TO 601 
DO 500 J =2, NLAST 
DO 501 1=1, NDMAX 

Y. I) =YY ( I , J ) 

CONTINUE 

CALL INFOPLT ( IC, NDMAX, X, 1 , Y, 1 , ZMIN, ZMAX, YMIN, YMAX, 1.0, 
1 1 , 1 H , 1 , 1 H , 0, SXC..SYC, 0.75,0. 75) ' 

CONTINUE 

CONTINUE 

PLOT LAST CURVE 

DO 60 1=1, NDMAX 
Y ( I ) =YY ( I , NG) 

CONTINUE 

CALL INFOPLT ( 1 , NDMAX, X, 1 , Y, 1 , ZMIN, ZMAX, YMIN, YMAX, 1 . 0, 

1 1,1H , 1 , 1 H , 0,SXC,SYC, 0.75, 0.75) 

CONTINUE 
CALL NFRAME 
RETURN 
CONTINUE 
GET DATA TO PLOT 
DO 701 1=1, NDMAX 
X ( I ) =FDPTDAT ( I , 1) 

DO 702 J = 1 , NG 
NN=JJJ— 4+ ( j-1) *5 
YY ( I , J ) =FDF'TDAT < I , NN) 

CONTINUE 

FDF'TDAT IS STORED IN THE FORM: 

Z , FLOWRATE , DENS I TY , PRESSURE , TEMPERATURE , Z , FLOWRATE , . . . 
CONTINUE 

FIND MAX AND MIN FOR Y VALUES 
YMIN=100. 

YMAX=0 . 0 

DO 703 1=1, NDMAX 

DO 704 J= 1 , NG 

I F ( YY ( I , J ) .LT. YMIN) YMIN=YY(I,J> 

IF ( YY ( I , J ) .GT. YMAX) YMAX=YY(I,J) 

CONTINUE 
CONTINUE 
MAX=10. *YMAX+1 . 0 
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705 

C 

c 


800 

C 


801 

C 


802 

C 


820 

C 
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Y, 1,ZMIN, ZMAX, YMIN, YMAX, 
SCALED , 

. 0 . 


504 


503 

5031 

C 


505 


YMAX=FLOAT (MAX) /10. 

MIN=10. *YM I N- 1 . 0 
YMIN-FLOAT (MIN) / 10. 

ZMIN=0. 0 
ZMAX=1.0 

PLOT FIRST CURVEC 
DO 705 1=1, NDMAX 
NN= 1 

Y ( I ) =YY ( I , NN ) 

PLOT JJJ V3 Z 
IF ( JJ J . GT. 6) GO TO 800 
JJJ=6, PLOT FLOWRATE VS Z 

CAi_L INFOPLTdC, NDMAX, X, 1, Y, 1 , ZM I N , ZMAX , YM I N, YMAX , 1.0. 

1 25, 25HZ , AXIAL DISTANCE, 

2 15, 15HFL0WRATE CM/SEC 

3 1 0 . , 5 . ,0.75,0.75) 

GO TO 820 
CONTINUE 

IF (JJJ .GT. 7) GO TO SOI 
JJ J=7, PLOT DENSITY VS Z 

CALL INFOF'LT ( IC, NDMAX, X, 1 , Y, 1 , ZMIN, ZMAX, YMIN, YMAX, 1 . 0, 

1 25, 25HZ , AXIAL DISTANCE, SCALED ,' 

2 16, 16HDENSITY, KG/M**3 ,0, 

3 10. , 5 . ,0.75,0.75) 

GO TO 820 

CONT I NUE 

IF (JJJ .GT. 8) GO TO 802 
J J J=8, PLOT PRESSURE VS Z 

CALL INFOPLTdC, NDMAX, X, 1,Y, 1 , ZMIN, ZMAX, YMIN, YMAX , 1.0, 

1 25, 25HZ , AXIAL DISTANCE, SCALED , 

2 14, 14HPRESSURE, TORR , 0 , 

3 10. ,5. ,0.75,0.75 ) 

GO TO 820 
CONTINUE 

J J J=9 , PLOT TEMPERATURE VS Z 

CALL INFOPLTdC, NDMAX, X, 1,Y, 1 , ZMIN, ZMAX, YMIN, YMAX , 1.0, 

1 25 , 25HZ , AXIAL DISTANCE, SCALED , 

2 IS, 1 8HTEMPERATURE , DEG C , 0 , 

3 10. ,5. ,0.75,0.75 ) 

GO TO 901 
PLOT REST OF THE CURVES 
IF (NG .EQ. 2) GO TO 5031 
DO 503 J=2, NLAST 
DO 504 1=1, NDMAX 
NN=J 

Y ( I ) =YY ( I , NN) 

CONTINUE 

CALL INFOF'LT (IC, NDMAX, X, 1 ,Y, 1, ZMIN, ZMAX , YMIN, YMAX , 1. 

1 1,1H , 1 , 1 H ,0, 10. ,5. ,0.75,0.75 ) 

CONTINUE 

CONTINUE 
PLOT LAST CURVE 
DO 505 1=1, NDMAX 
Y ( I ) =YY (I , NG) 

CONTINUE 


: 1 0 . , 5 • , C 

>.75, 

0.7 

CONTINUE 



IF (NLAST 

.EQ. 

0) 


0 , 
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CALL I NFOF’LT ( 1 , NDMAX , X , 1,Y, 1, ZMIN, ZMAX , YMIN, YMAX , 1.0, 
1 1,1H , 1,1H ,0, 10. ,5. ,0.75,0.75 ) 

90 1 CONT I NUE 

CALL NFRAME 

RETURN 

END 


SUBROUTINE F'FLOW 

C SUBROUTINE TO CALCULATE THE PARAMETERS FOR SUBROUTINE FLOW 

C PARAMETERS STORE IN COMMON BLK10 

COMMON /BLK7/ ABC , COO, CO, 0MEG1 , P, R1 , R2, TM, XNRHO 
COMMON/ BL.K8 / Z OL , ZE, NG, TO, RAD, A, PIN, W 
COMMON/'BLK 1 0/CF1 , CF2, CF4, QFO, R3TAR, ZL, L, SFi , SF2 , AR , AAO , BBO 
C0MM0N/BLK23/W0, ETAO, PTO, FRAC 
C OMMON / BLK29 / Z Z Z , TZ Z , P Z Z , ET A Z Z , W Z Z 
REAL L 
C 

C AR IS BEAM DIAMETER RADIUS IN CM 

S- 1=1000. / 296 . 

C SFI =1 / <M*1 . E-03) , WHERE M IS THE MOLECULAR WEIGHT 

C C/s (J/K KG) =SF1*CV (J/K MOL) 

SF 2= ( 1 . 01325E5) /760. 

C SF2 EQUALS 133.3 F'A/TORR 

C SF 1 , SF2 ARE SCALE FACTORS FOR THE CORRECT UNITS OF 

C F -PRESSURE <N/M**2> , SFI (MOLE/KG) 

C F'T=PRESSURE (TORR), SF2 (N/M**2) /TORR 

C TEST TO SEE IF WO <M/SEC) IS NEAR MAX VALUE OF 5. 30*SQRT (TO) 

TEST=5 . 3*SQRT (TO) 

TEST 1 =ABS ( 100. * (WO-TEST) /TEST) 

IF (TEST 1 .LT. 10.0) WRITE<6,778) 

77S FORMAT ( IX, 52HWARNING WO VALUE IS WITHIN 10 PERCENT OF ITS MAXIMUM 
1 , /, IX, 35HALL0WABLE VALUE OF 5.3*SQRT(T0) ) 

C ZL=LIGHT SOURCE LENGTH IN CM 

C ETA=DENSI TY <KG/M**3> 

C L=TUBE LENGTH IN CM 

C RSTAR=GAS CONSTANT (JOULE/KG DEG K) 

C T=TEMF'ERATURE (DEG K) 

C CV=SPECIFIC HEAT AT CONSTANT VOLUME 

C RSTAR (J/K KG) = SF1*R<J/K MOL) 

C WO, WL=FLOW VELOCITY (M/SEC) (SUBSCRIPTS 0,L FOR START, END) 

PO=SF2*PTO 
RST AR=83 14. 3/296.0 
CF1=ETA0*W0 
CF2=CF1*W0+P0 

QFO=FRAC* ( 1 . 40E3) *C0 *2. / (A*l . E-2) 

C QFO IN W/M**3 

C INPUT (W/M**2) DISTRIBUTED OVER VOLUME (2 (PI ) A) /PI < A**2) IN M2/M3 

WZZ=WO 
RETURN 
END 


SUBROUT I NE FLOW ( Z , T, PTORR , WSS , ETA , DWDZ ) 
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COMMON /BLK8/Z0L , ZE, NG, TO, RAD, A, PIN, W 
C SUBROUTINE TO CALCULATE T,F',WSS,ETA AS FUNCTION OF Z 

COMMON/BLKIO/CF 1 , CF2, CF4 , QFO, RSTAR, ZL, L, SF1 , SF2, AR, AAO, BBO 
COMMON/BLK23/WO , ETAO, PTC, FRAC 
C 

COMMON /BLK29/ ZZZ, TZZ, PZZ , ETAZZ , WZZ 
REAL L 
I C0UNT=0 
Q=QFO*Z/ 100. 

C ENERGY INPUT TERM DETERMINED BY QFO 

C Z IS IN CM AND QFO IS IN W/M**3 

IFCZ.GT. ZL) Q=QF0*ZL/100. 

X X X =BB0* ( T 0-300 . ) 

CALL ETO i XXX , EEC) 

WSS=WZZ 

50 CONTINUE 

T = ( (CF2/CF1 ) *WSS-WSS*WSS ) /RSTAR 
DTDW= ( (CF2/CF1 ) -2. *WSS) /RSTAR 
XXX--BBO* (T-300. ) 

CALL ETO (XXX,, EE1 ) 

FIRSTAR* (T-TO) + ( AAO/BBO ) *SF1* (EEl-EEO) + 

1 . 5* (WSS*WSS~WO*WO> -Q/CFl 
FP= ( RST AR+SF 1 * AAO*EE 1 ) *DTDW+WSS 
W1=WSS-F/FF 

ERR=ABS (100. * (Wl-WSS) /WSS) 

IF (ERR -LT. .25) GO TO 100 

IF ( I COUNT . GT. 95) WRITE <6, 357) ICQUNT, Z , WZZ , WSS, CF1 , CF2, TO, Q, FM 
357 FORMAT ( IX, 15, IX, (3E16. 7, /) ) 

WSS=W1 

IC0UNT=IC0UNT+1 
IF (I COUNT .GT. 100) STOP 4444 
GO TO 50 
100 CONTINUE 
WSS=W1 

T=( (CF2/CF1 ) *WSS-WS5*WSS ) /RSTAR 
TC=T— 273 
ETA=CF 1 /WSS 
PNM2=ETA*RSTAR*T 
PT0RR=PNM2/SF2 

112 = 2 

TZ Z=TC 
PZZ=F'TORR 
ET AZ Z=ET A 
WZZ=WSS 

XX X=BBO* (T-300. ) 

CALL ETO (XXX, COO) ^ 

CV=CVO*AAO 
CVS=SF 1 *CV 
XNUM=RST AR*QFO 

XDEN= (CF2-2. *CF 1*WSS) * (CVS+RSTAR) +RSTAR*CF1*WSS 
DWDZ=XNUM/XDEN 
RETURN 
END 
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SUBROUTINE ARREN ( TEMP ) 

SUBROUTINE FOR ARRENHIUS EXPRESSION OF RATE COEFFICIENTS 
BASIC ASSUMPTIONS 

FOR Q I TERMS Q I =0 I 0*EXP ( -BETA* ( TEMP-TO) ) 

TREAT KI TERMS LIKE Cl TERMS 

C0MM0N/BLK2/K 1 , K2,K3, K4 , K5 , K6 , K7, KB, Cl , C2, C3,C4, C5 
C0MM0N/BLK27/Q 1 , Q2,Q3,Q4, 05 

COMMON / BLK3 / B , B2 , B3 , C , AOO , BOO , EPSNU, OMEGA, C6 
COMMON /BL.K1 1 /KK1 , KK2 , KK3 , KK4 , KK5 , KK6 , KK7 , KK3 
COMMON /BLK1 2 /QQt. , 002,003, 004,005 
CCMM0N/BLK13/CC 1 , CC2, CC3, CC4, CC5, CC6 
REAL Kl , K 2 , K3, K4,K5, K6, R7, K8, KK1 , KK2, KK3, KK4, KK5, KK6, KK7, KK8 
C 

C REFERENCE J.S. COHEN AND O.P. JUDD 

C J.APF'L. RHYS., VOL 55, NO. 7, APRIL 1984 

C COEFFICIENTS MODIFIED TO ACHIEVE SPECIFIC VALUES AT TEMFERATURE 

C -OF 276 DEGREES K. 

C 

BETA=4 . 4E-3 
SF 1=1.0 

X X X=~BET A* ( TEMP-300 ) 

CALL ETO (XXX, YYY) 

SF2=YYY 
K 1 =KK 1 *SF 1 
K2=KK2*SF 1 
K3=KK3*SF 1 
K4=KK4*SF1 
K5=KK5*SF1 
K6=KK6*SF1 
K7=KK7*SF 1 
K8=KK8*SF 1 
Cl =CC1 *SF 1 

C2=CC2*EXF < 1360. 00/TEMP) 

C3=CC3*SF 1 

XYZ=-29. 437-5. 844*AL0G10 (TEMP/300. )+2. 163* (AL0G10 (TEMP/300. > > **2 
C4=CC4**XYZ 

C5=CC5*EXP ( 1310. 000 /TEMP) 

C6=CC6*SF1 
Q1=0Q1*SF 1 

02=002*EXP( (-4. 4E-3)* (TEMP-300. ) ) 

Q3=QQ3*SF 1 
04=QQ4*SF 1 
Q5=QQ.5*SF 1 
RETURN 
END 


SUBROUTINE ETO(X,Y) 

C NEGATIVE EXPONENTIAL FUNCTION 

IF ( X .LT. -670.) GO TO 100 

ORIGINAL page is 

OF POOR QUALITY 
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100 


Y=EXP ( X ) 
RETURN 
Y=0 . 
RETURN 
END 


FUNCTION CHS I 1 (Z) 

Z 1 BAR IS CUTOFF POINT OF ILLUMINATION 
ABARO IS FRONT CUTOFF POINT OF ILLUMINATION 
CHS II IS A CONSTANT 
IMPLICIT REALMS ( A-H , K , L , 0- Z ) 

CGWM0N/BLK4/CHSI 10, CHS I 20, ABARO, Z 1 BAR , LC , XKO 
COMMON/BLK8/ZOL, ZE, NG, TO, RAD, A, PIN, W 
REAL LC 

IF (Z. LT. ABARO) GO TO 100 
IF (Z.LT. Z1BAR) GO TO 200 
Z GREATER THAN Z 1 BAR 
00 CHS 1 1 =0 . 0 

CHS II HAS UNITS OF SEC'-l 
RETURN 
CONTINUE 
CHS 1 1 =CHS I 1 0 
RETURN 
END 


FUNCTION CHS I 2 ( Z ) 

ABARO IS FRONT CUTOFF POINT FOR ILLUMINATION 
Z1BAR IS CUTOFF POINT FOR ILLUMINATION 
CHS I 2 IS A CONSTANT 
IMPLICIT REALMS (A-H, K,L,0-Z> 

COMMON/BLK4/CHSI 10, CHS 1 20, ABARO, Z 1 BAR, LC, XKO 
COMMON / BLK8 / ZOL , Z E , NG , TO , RAD, A, PIN, W 
REAL LC 

IF (Z.LT. ABARO) GO TO 100 
IF (Z.LT. Z1BAR) GO TO 200 
C Z GREATER THAN Z 1 BAR 

CC XXX=- (Z-Z1BAR) /ZE 

CC CALL ETO ( XXX , YYY) 

CC CHS 1 2=CHS 1 20*YYY 

C CHSI2 HAS UNITS OF SEC-"-! 

RETURN 

100 CHS I 2=0.0 

RETURN 

200 CONTINUE 

CHS I 2=CHS I 20 

RETURN 

END 
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SUBROUTINE COEFFS 

THIS SUBROUTINE DEFINES THE COEFFICIENTS IN THE DIFFERENTIAL 
EQUATIONS TO BE SOLVED. 

IMPLICIT REAL*8 ( A-H, K, L, O-Z ) 

COMMON/BL.K22/AD, V 1 , V2 , GG 

C0MM0N/BLK3/B, B2, B3, C, AGO, BOO, EF'SNU, OMEGA, C6 
COMMON / BLK4/CHS 110, CHS I 20 , ABARO , Z 1 BAR , LC , XKO 
COMMON /BLK7/ ABC, COO, CO, 0MEG1 , P, R1 , R2, TM, XNRHO 
COMMON / BLK8 / ZOL, ZE, NG, TO, RAD, A, PIN, W 
C3MM0N/BLK1 l /KK1 , KK2, KK3, KK4, KK5 , KK& , KK7 , KK8 
COMMON /BLK12/QQ1 , QQ2, QQ3, QQ4 , QQ5 
COMMON/ BLK 1 3/CC 1 , CC2 , CC3, CC4, CC5, CC6 
COMMON /BLK28/T AU2 , T AU4 , T AU5 , T AU6 , T AU3, S IG 
C0MM0N/BLK299/TTT2 , TTT3 , TTT4, TTT5, TTT6 
REAL KK 1 , KK2 , KK3 , KK4 , KK5 , KK6 , KK7 , KK8 , LC 

COEFFICIENTS IN THE DIFFERENTIAL. EQUATIONS 

0MEGA=0MEG1 

OMEGA=OMEG 1 IS FLOW RATE IN CM/SEC 
OMEGA IS FLOWRATE AT RADIUS R=RAD 
SEE SUBROUTINE VELOC 

CALL VELOC (0MEG1 , RAD, OMEGA, A) 

ABAR0=0 . 0 

ABARO IS START OF ILLUMINATION 

Z 1 BAR=2*Z0L IS THE POINT ON AXIS WHERE ILLUMINATION BEGINS TO DIMINISH 
CHS I 10= <3. 04E-3) *CO*XNRHO 
C COO IS THE CONCENTRATION IN SOLAR CONSTANTS 

CHSI20= <3. 38E-2) *C0/2. 5 
Z1BAR=2*Z0L 
EPSNU= 1 . 5E-19 
C WATTS*SEC 

AOO=2. 0E17 
BOO=0. 443 
TAU2=TTT2 
TAU3=TTT3 
TAU4=TTT4 
TAU5=TTT5 
TAU6=TTT6 
C 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 


c 

(CHSI1) 

HV + 

RI — 

> R + I* 


c 

(CHSI2) 

HV + 

12 — 

> I* + I 

(517. OF THE TIME 

c 

p 

(CHSI2) 

HV + 

12 — 

> I + I 

(497. OF THE TIME) 

u. 

c 

( KK 1 ) 

I* + 

R — 

> RI 


c 

<KK2> 

I + 

R — > 

RI 


c 

(KK3) 

R + 

R — > 

R2 


c 

(KK4) 

R + 

RI — 

> R2 + I 


c 

<KK5) 

R + 

12 — 

> RI + I 


c 

(KK6) 

R + 

RI — 

> R2 + I* 


c 

(KKS) 

I *++RS I — 

>>II2++RR 



c 
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c 

< OQ 1 ) 

I* + R I — > I 

+ RI 

c 

(QQ2) 

I* + 12 — > I 

+ 12 

c 

(QQ3) 

I* + R — > I 

+ R 

c 

(QQ4) 

I* + R2 — > I 

+ R2 

c 

(QQ5) 

I* + I — > I 

+ I 

c 




cc 

(CC1) 

I* + I + RI 

— > 12 + RI 

cc 

(CC2) 

I+I+RI- 

-> 12 + RI 

cc 

(CC3) 

I* + I + 12 

— > 12 + 12 

cc 

(CC4 ) 

I + I + 12 - 

-> 12 + 12 

cc 

(CCS) 

I + I + R2 - 

-> 12 + R2 

cc 

(CC6) 

I + I * + R2 

— > 12 + R2 

c 




c 




c 

WALL REACTIONS 


c 




c 

(VI ) 

I + I + WALL — 

■> 12 + WALL 

c 

<V2) 

R2 + I + WALL — 

-> R + RI + WALL 

c 




cccccccc 

ccccccccccccccccccccccccccccccccccccccccc 


c 

c 

c 

C WILSON ET AL. (1984) AND STOCK ET AL. (1886) 

C KK1 * 5.6E-13 

C COHEN AND JUDD (1984) 

C KK1 = . 9E-13 

C V I NOKKURQ V AND ZALESSKKII (1979) 

KK1 = l.E-14 
C 
C 

C WILSON ET AL. (1984) AND STOCK ET AL. (1986) 

C KK2 = 2.3E-11 

C BREDERLOW ET AL. (1983) 

KK2 = . 8E- 1 1 
C COHEN AND JUDD (1984) 

C KK2 = . 7E-1 1 

C 
C 
C 

C WILSON ET AL. (1984) AND STOCK ET AL. (1986) 

C KK3 = 2.6E-12 

C BREDERLOW ET AL. (1983) AND COHEN AND JUDD (1984) 

KK3 * 2.0E-12 
C 
C 

C WILSON ET AL. (1984)., STOCK ET AL. (1936), AND COHEN AND JUDD (1984) 
KK4 = 3.E-16 
C 
C 

C COHEN AND JUDD (1984) 

KK5 = l.OE-U 
C 

C WILSON ET AL. (1984) AND STOCK ET AL. (1986) 

C KK6 = 3.2E-17 
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C BREDERLOW ET AL. (1983) AND COHEN AND JUDD (1984) 

KK6 = 0. 

C 

C 

C COHEN AND JUDD (1984) 

KK7 = 3.0E-19 
C 
C 
C 

C COHEN AND JUDD (1984) 

KK3 = 1.6E-23 
C 
C 

C WILSON ET AL. (1984) AND STOCK ET AL. (1986) 

C QC1 = 2.0E-16 

C COHEN AND JUDD (1984) 

C QQ 1 = 1.7E-16 

C BREDERLOW ET AL. (1983) 

QQ 1 = 1.7E-17 
C 
C 
C 

C COHEN AND JUDD (1984) 

C QQ2 = 3.80E-1 l*EXP<-4. 4E-3* (T-300. ) ) 

QQ2=3. 8E- 1 1 

C SUBROUTINE ARREN ADDS TEMPERATURE EFFECT 

C BREDERLOW ET AL. (1983) 

C QQ2 = 3.0E-11 

C WILSON ET AL. (1984) AND STOCK ET AL. (1986) 

C QQ2 = 1.9E-11 

C 
C 
C 

C COHEN AND JUDD (1984) 

QQ3 = 3.7E-18 
C 
C 

C COHEN AND JUDD (1984) 

QQ4 = 4.7E-16 
C 
C 

C COHEN AND JUDD (1984) 

QQ5 = 1.6E-14 
C 
C 
C 

C COHEN AND JUDD (1984) (A DIFFERENT REACTION, PERHAPS) 
CC CC1 = l.E-32 

CC WILSON ET AL. (1984) AND STOCK ET AL. (1986) 

CC CC1 = 3.2E-33 

CC HOHLA AND KOMPA (1976) 

CC1 = 1.6E-33 
CC 
CC 

C COHEN AND JUDD (1984) 

C CC2 = 5. 7E-33*EXP ( 1360. /T> 
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CC2 = 5.7E-33 

C SUBROUTINE ARREN ADDS TEMPERATURE EFFECT 

CC HOHLA AND KOMPA (1976) 

CC CC2 = 2. 1E-33*EXP ( 1600. /T) 

C BREDERLOW ET AL. (19B3) 

CC CC2 = 3.8E-31 

C WILSON ET AL. (1984) AND STOCK ET AL. (1986) 

CC CC2 = 8.5E-32 

C 

c 

C WILSGN ET AL. (1984) AND STOCK ET AL. (1986) 

CC CC3 = 8.E-32 

C BREDERLOW ET AL. (1983) AND COHEN AND JUDD (1984) 
CC3 = 0. 

C 

C 

C WILSON ET AL. (1984) AND STOCK ET AL. (1986) 

CC C24 = 3.8E-30 

C COHEN AND JUDD (1984) 

C CC4 = 10. 0** (-29. 437-5. B44*AL0G10 (T/300. ) 

CC4 = 10.0 

C SUBROUTINE ARREN ADDS TEMPERATURE EFFECT 

C BREDERLOW ET AL. (1933) 

CC CC4 = 2.9E-30 

C 

C 

C COHEN AND JUDD (1984) 

C CCS = 8. E-33*EXP ( 1310. /T) 

CCS = 8.0E-33 

C SUBROUTINE ARREN ADDS TEMPERATURE EFFECT 

C 

C 

C BREDERLOW ET AL. (1983) AND COHEN AND JUDD (1984) 
CC6 = 0.0 
C 
C 

C VI = 0.0 

C 

C 

C V2 = 0.0 

C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c 

c 

C KK 1 = . 903E- 1 3 

C KK2=80.E— 12 

C KK3=. 65E-12 

C KK4=1 . 000E-16 

C KKS=3.0B9E?1 1 

C KK7=. 1517E-1S 

C KK8=1.6E-23 

C 

C QQ1=. 476E-16 

C QQ2=1 . 9E-1 1 

C QQ3=. 1235E-17 

C QQ4= 1 . 57E- 1 6 
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C 005= . 53E- 1 4 

C 

CC CC 1=1. 053E-33 

CC CC2=45. OE-32 

CC CC3=. 4447E-31 

CC CC4=4. 94E-30 

CC CC5=3. 6E-3 1 

CC CC6= 1.8E-32 

C 

C V 1 = 1.0E-12 

C V2= l.OE-11 

C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c 

c 

GG = 2* ( . 18/LC) **2 
C=3. 0E10 

CC B=F’* <3. 5E16) 

B= (9. 66E18) *P/TO 

B2=B*B 

B3=B2*B 

RETURN 

END 


SUBROUTINE VELOC <OMEGl , RAD, OMEGA, A) 

C VELOCITY PROFILE IS ASSUMED TO BE PARABOLIC 

C MAX VELOCITY AT RAD-0.0 

C ZERO VELOCITY AT RAD-A 

C A IS RADIUS OF TUBE 

C OMEG1 IS MAXIMUM VELOCITY ALONG CENTERLINE 

C 

C CALCULATE VELOCITY OMEGA AT R=RAD 

C 0 .LE. RAD .LE. A 

C 

C TYPE OF FLOW 

OMEGA= <OMEGl / (A* A) ) * (RAD-A) **2 

RETURN 

END 


SUBROUTINE PRINT ( Z , TEMP) 

COMMON/BLK2/K1 , K2,K3, K4, K5, K6, K7, K8, Cl , C2, C3, C4, C5 
C0MM0N/BLK27/Q 1 ,02,03 , 04 , Q5 

C0MM0N/BLK3/B, B2, B3, C, AGO, BOO, EPSNU, OMEGA, C6 
COMMON /BLK4 /CHS 1 10, CHS 1 20, ABARO, Z 1BAR, LC, XKO 
COMMON /BLK7 / ABC, COO, CO, 0MEG1 , P, R1 , R2, TM, XNRHO 
COMMON/BLK22/AD, VI , V2, GG 
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c 

10 

1 1 
c 


1 00 

101 

102 

103 

104 

105 


C0MM0N/BLK28/TAU2, TAU4, TAU5, TAU6, TAU3, SIG 
COMMON / BLK299 / TTT2 , TTT3, TTT4, TTT5, TTT6 
REAL K 1 , K2 , K3 , K4 , K5 , K6 , K7 , KB , L.C 

IF (Z . EQ. 0.0) WRITE (6, 10) 

FORMAT ( IX, 20HC0EFFICIENTS AT 2=0 ) 

IF ( Z .GE. 1.0) WRITE (6, 11) 

FORMAT (IX, 20HC0EFF I C I ENTS AT Z=L ) 

WRITE OUT COEFFICIENTS 
WR I TE ( 6 , 1 00 ) K 1 ,K7, Q1 , C5, TTT2 
WRITE ( 6 , 101) K2, K8, Q2, 01 , TTT3 
WRITE (6, 102) K3, Cl , Q3, 02, TTT4 
WRITE < 6 , 103) K4 , C2 , Q4 , S I G , TTT5 
WR I TE (6,1 04 ) K5 , C3 , Q5 , TEMP , TTT 6 
WR I TE ( 6 , 1 05 ) K6 , C4,C6, LC 

FORMAT <T5,5HKK1 » , El 5. 7, T30, 5HKK7 = , E15. 7, T60, 5HQQ1 = 
1 T85,5HCC5 = , E15.7 , T103,7HTTT2 = ,E15.7 ) 

FORMAT (T5,5HKK2 = , E15. 7 , T30, 5HKK3 = , El 5. 7, T60, 5HQQ2 = 
1 T85, 5HV1 = , E15.7 , T103,7HTTT3 = ,E15.7) 

FORMAT (T5, 5HKK3 = , El 5. 7, T30, 5HCC1 = , El 5. 7, T60, 5HQQ3 = 
1 T85, 5HV2 = , E15.7 , T103,7HTTT4 = ,E15.7 ) 

FORMAT (T5, 5HKK4 = , El 5. 7, T30, 5HCC2 = , E15. 7, T60, 5HQ04 = 
1 785, 8HSIGMA = , E15.7 , T103,7HTTT5 = ,E15.7 ) 

FORMAT (T5, 5HKK5 = , El 5. 7, T30, 5HCC3 = , E15. 7, T60, 5HQQ5 = 
1 T85, 7HTEMP = ,F9.2 , T103,7HTTT6 = ,E15.7 ) 

FORMAT <T5, 5HKK6 = , El 5. 7, T30, 5HCC4 = , E15. 7, T60, 5HCC6 = 
1 T85, 4HL = , F9 . 2 ) 

RETURN 

END 


SUBROUTINE FUN(N,S,Y,F) 

THIS SUBROUTINE DEFINES THE RIGHT HAND SIDE 
OF THE DIFFERENTIAL EQUATIONS FOR THE CHEMICAL KINETICS 
IMPLICIT REAL*8 ( A-H, K, L, 0-Z ) 

DIMENSION Y (7) ,F(7) 

COMMON/BLK1/X7, POWER 
EXTERNAL CHS II, CHS 1 2 

COMMON / BLK2 / K 1 , K2 , K3 , K4 , K5 , K6 , K7 , K8 , C 1 , C2 , C3 , C4 , C5 
COMMON/ BLK27 / Q 1 , Q2, Q3, Q4, Q5 

COMMON / BLK3 / B, B2, B3, C, AOO, BOO, EPSNU, OMEGA, C6 
COMMON/BLK4/CHSI 10, CHS 1 20, ABARO, Z 1 BAR, LC, XKO 
COMMON/BLK7/ABC, COO, CO, 0MEG1 , F', R1 , R2, TM, XNRHO 
C0MM0N/BLK22/AD, VI , V2, GG 

COMMON/BL.K28/TAU2, TAU4, TAU5, TAU6, TAU3, SIG 
C0MM0N/BLK299/TTT2 , TTT3, TTT4, TTT5, TTT6 
REAL K 1 , K2 , K3 , K4 , K5 , K6 , K7 , K8 , LC 

QY=QUANTUM YIELD 
Q Y =1.0 


F ( I ) , 1=1 , 6 ARE RATES OF CHANGES FOR THE CONCENTRATIONS 
F ( 1 ) =DCRI ] /DZ F (2) =DCR1/DZ 

F(3)=DCR2I/DZ F (4) =DC 123 /DZ 


E15.7 , 
E15.7 , 
El 5. 7 , 
E15.7 , 
E15.7 , 
E15.7 , 
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F (6) =DC 1 1 /DZ 


F(5)=DCI*:/DZ 

0 CONTINUE 
Z=LC*S 

S IS SCALED LENGTH VARIABLE 
0 . LE , S . LE . 1 

FLOW CALCULATES TEMPERATURE,, PRESSURE, FLOWRATE 
AND DENSITY RESPONSE AS FUNCTION OF Z 

Z IS DISTANCE IN CM 

CALCULATE GAS PARAMETERS AS FUNCTION OF Z 
CALL FLOW ( Z, TEMP, PRESS, FLOWR, DENSITY, DWDZ ) 

CALL FLOW ( Z , TEMP , PRESS , FLOWR , DENSITY , DWDZ ) 

WRITE (6, 357) Z, TEMP, PRESS, FLOWR, DENSITY, DWDZ 
57 FORMAT ( 1 X , 3E 16.7, / , 3E 16.7 ) 

TEMP IS TEMPERATURE DEG K 
PRESS IS PRESSURE IN TORR 
FLOWR IS FLOWRATE IN M/SEC 
DENSITY IS GAS DENSITY IN KG/M**3 
OMEGA=FLOWR* 100. 

OMEGA IS FLOW RATE IN CM/SEC 

CALCULATE COEFFICIENTS AS FUNCTION OF TEMP AND Z 
CALL ARREN (TEMP) 

TAU5- (TTT5) *PRES8 
TAIJA— TAU5 
C 
C 

C CONSTANTS COME VIA COMMON BLKS 2 AND 3 

C K’S IN CM**3/SEC 

C C ! S IN CM**A/SEC 

C Q’S IN CM**3/SEC 

C RBWRROZSYIBJ) MBCM**2 

X7STAR=Y(7)*B+X8 
DIF=Y(5)-.5*Y(A) 

CALL SIGMA (SIG2) 

SIG=SIG2 

F ( 1 ) =K1*B*Y (2) *Y (5) +K2*B*Y (2) *Y (6) -CHS 1 1 (Z) *Y ( 1 ) -K4*B*Y ( 1 ) *Y (2) 

1 +K5*B*Y(2)*Y<4> -K7*B*Y (5) *Y ( 1 ) -KA*B*Y (2) *Y ( 1 ) +V2*B*Y (3) *Y (6) 

2 -K8*B*Y<A)*Y<1) -Y ( 1 ) *DWDZ 

F (2) =CHSI 1 (Z) *Y ( 1 > -K1*B*Y <2) #Y (5) -K2*B*Y ( 2 ) *Y (A) -2*K3*B*Y (2) *Y (2) 

1 -K4*B*Y ( 1 ) *Y (2) -KA*B*Y ( 1 ) *Y (2) -K5*B*Y (2) *Y <4> +V2*B*Y (3) *Y (6) 

2 +K7*B*Y (5) *Y ( 1 ) +K8*B*Y (6) *Y ( 1 ) -Y (2) *DWDZ-Y (2) /TAU2 

F ( 3 ) =K3*B*Y ( 2 > *Y ( 2 > +KA*B*Y < 1 > *Y ( 2 > +K4*B*Y ( 1 > *Y ( 2 ) - V2*B*Y (3) *Y (A) - 
1 Y (3) *DWDZ -Y(3)/TAU3 

A1=C1*B2*Y ( I ) *Y (5) *Y (A) +C2*B2*Y ( 1 ) *Y ( A) *Y (A) +C3*B2*Y <4) *Y <5) *Y <A> 
A2=C4*B2*Y (4) *Y (A) *Y (A) -CHS 1 2 (Z) *Y <4) +K7*B*Y(5)*Y ( 1 ) 

1 -K5*B*Y(2)*Y(4) +V1*B*Y<A)*Y(A> +C5*B2*Y (A) *Y (A) *Y (3) 

F f 4) =A1+A2+K8*B*Y (A) *Y ( 1 ) +CA*B2*Y (A) *Y (5) *Y (3> -Y (4) *DWDZ-Y (4) /TAU4 
A3=QY*CHSI 1 <Z ) *Y < 1 ) +0. 51*CHSI2 < Z> *Y <4> -K1*B*Y (2) *Y <5) 

A4=-C 1 *B2*Y ( 1 > *Y (5) *Y (A) -C3*B2*Y (4) *Y <5) *Y (A) -Q1*B*Y < 1 ) *Y (5) 
A5=-Q2*B*Y (4) *Y (5) -C*SIG*X7STAR*DIF +KA*B*Y <2> *Y ( 1 > 

F (5) =A3+A4+A5-Q3*B*Y (5) *Y (2) -Q4*B*Y (5) *Y (3) -Q5*B*Y (5)*Y (A) 

1 -K7*B*Y (5) *Y < 1 ) -CA*B2*Y (A) *Y <5> *Y (3) -Y (5) *DWDZ-Y <5) /TAU5 
AA= 1 . 49*CHSI2 <Z ) *Y <4> +Q1*B*Y ( 1 ) *Y (5) +Q2*B*Y (4) *Y (5> 

1 -2*C5*B2*Y (A) *Y (A) *Y (3) -K8*B*Y < A) *Y < 1 ) 

A7=C*SIG*X7STAR*DIF -C1*B2*Y ( 1 ) *Y (5) *Y (A) 
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A8=-2*C2*B2*Y ( 1 ) *Y (6) *Y (6) -C3*B2*Y (4) *Y <5) *Y (6) 

A9=-2*C4*B2*Y (4) *Y <6> *Y (6> -K2*B*Y (2) *Y (6) +K4*B*Y < 1 ) *Y (2) 
A10=Q3*B*Y (5) *Y (2) +Q4*B*Y (5) *Y <3) +Q5*B*Y (5) *Y (6) 

1 +K5*B*Y (2) *Y (4) -V2*B*Y < 3) *Y ( 6 > -2*V1*B*Y (6) *Y (6) 

F (6) =A6+A7+A8+A?+A10-C6*B2*Y (6) *Y (5) *Y (3) -Y (6) *DWDZ-Y (6) /TAU6 
SCALED EQUATIONS IN THE Z-DIRECTION 
DO 10 1=1,6 
F(I)=LC*F(I> /OMEGA 
F (7) =LC*Y (7) *DIF*B*SIG 
RETURN 
END 


SUBROUTINE SIGMA (SIG) 

THIS SUBROUTINE DEFINES THE CROSS SECTION SIGMA 
COMMON/ BLK3/B , B2 , B3, C, AOO, BOO, EF'SNU, OMEGA, C6 
COMMON /BLK7/ABC, COO, CO, 0MEG1 , P, R1 , R2, TM, XNRHO 
REAL NU , NUS , NUO , NU 1 , NU2 , NU3 , NU4 , NU5 
P I =3 . 14159 
F I S=P I *P I 
NU=C/ 1 . 315246E-4 
NUS=NU*NU 
P I SNUS=P I S*NUS*4 . 

G=0 

cs=c*c 

NUO=NU 

NU1=NU0+. 141*C 
NU2=NU1+. 068*0 
NU3=NU0— . 427*C 
NU4=NU3- . 026*C 
NU5=NU4-. 068*C 
DELT A23=NU-NU5 
DELT A22=NU-NU4 
DELT A2 1 =NU-NU3 
DELT A34=NU— NUO 
DELT A33=NU— NU 1 
DELT A32=NU-NU2 
TEMP0=293 
TWALL=TEMPO 
T 1=TWALL 
A=5. 434 

A1=A*2. 4/7. 7*CS 
A2=A*3. 0/7. 7*CS 
A3=A*2. 3/7 . 7*CS 
A4=A*5. 0/7. 7*CS 
A5=A*2. 2/7. 7*CS 
A6=A*0. 6/7. 7*CS 
FUGTEMF-SQRT ( T 1 /300 . ) 

ALPHAM=1 . 8BE7*FUGTEMP 
DELD0P=2. 51E8*FUGTEMP 
DELNU=DELDOP+ALPHAM*P 

SIGMA23=A1 / (PISNUS*DELNU) / ( 1+ (2. *DELTA23/DELNU> **2> *5. /I 
SIGMA22=A2/ <F'ISNUS*DELNU> / ( 1 + (2. *DELTA22/DELNU> **2> *5. /I 
SIGMA21=A3/ ( P I SNUS*DELNU ) / ( 1+ (2. *DELTA21 /DELNU) **2> *5. /12. 
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SIGMA34=A4/ <PISNUS*DELNU> / ( 1+ (2. *DELTA34/DELNU> **2) *7./ 12. 

SIGMA33=A5/ (F'ISNUS*DELNU) / ( 1+ <2. *DELTA33/DELNU> **2) *7. / 12. 

S I GMA32=A&/ ( P I SNUS*DELNU ) / ( 1+ (2. *DELTA32/DELNU) **2> *7. / 12. 

S I GMAT=S I GMA23+S I GMA22+S I GMA2 1 +S I GMA34+S I GMA33+S I GMA32 

SIG=SIGMAT 

RETURN 

END 


SUBROUT I NE I NTEG ( I PR I NT , H ) 

THIS SUBROUTINE INTEGRATES THE SYSTEM OF DIFFERENTIAL EQUATIONS 
DEFINING THE CHEMICAL KINETICS OF N-C3F7I IODINE LASER 

USING A VARIABLE STEP SIZE 7TH ORDER RUNGE KUTT A-FEHLBERG METHOD. 
I MPL I C I T RE AL*8 ( A-H , K , L , 0- Z ) 

DIMENSION Y0<7) , X (7) , WK (84) 

COMMON /BLK1/X7, POWER 

COMMON/ BLK3/B, B2, B3, C, AOO, BOO , EPSNU , OMEGA , C£> 

COMMON / BLK4 / CHS 1 1 0 , CHS I 20 , ABARO , Z 1 BAR , LC , XKO 
COMMON /BLK7/ ABC , COO , CO , OMEG 1 , P , R 1 , R2, TM, XNRHO 
COMMON / BLK8 / ZOL, ZE, NG, TO, RAD, A, PIN, W 
COMMON /BLK22/ AD , VI , V2, GG 
CCMM0N/BLK23/W0’ ETAO, F'TO, FRAC 
COMMON /BLK29 / Z Z Z , T Z Z , P Z Z , ET A Z Z , W Z Z 
COMMON /BLK30/ DAT A (50,50) , NDMAX , FLRATE (8) , FDF'TDAT (50, 40) 

EXTERNAL FUN, CHSI 1 , CHSI2 
REAL LC 

INTEGRATE SYSTEM FROM S=0 TO S-1.0 USING RUNGE-KUTTA METHOD 

X (1)=RI 
X (2) =R 
X (3) =R2 
X (4 ) = 1 2 
X (5>=I* 

X (6)=I 
X (7) =RHO+ 

X8=RH0- 
X9= I *- . 5*1 

INITIALIZE CONSTANTS FOR FLOW EQUATIONS 

SEE COMMON BLK10 FOR THESE CONSTANTS — NEEDED FOR SUB FLOW 

CALL F'FLOW 

ND=0 

TEST FOR PRINT CONDITIONS 
IF ( IPRINT. EQ. 0) GO TO 229 
NG=NG+1 

NG IS THE NUMBER OF GRAPHS TO BE SAVED 
MAXIMUM NG=B 

FLOWRATE (NG) IS LABEL FOR DATA SAVED 

SEE ALSO PROGRAM PLOTD — WHICH CAN PLOT THE SAVED DATA 

FLRATE (NG)=OMEG A 
WRITE (8, 331 ) LC 
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331 

229 


C 

C 


C 


9 

C 


C 


C 

C 


191 


300 

lO 

C 


c 


222 

C 

C 

C 


FORMAT (IX, 1 1H LC = ,F10.2> 

CONTINUE 

CONTINUE 

N=7 

TQL=1 . OE-5 
PD=1 . 0 
MTH= 1 

wzz=wo 

H IS STEP SIZE IN SCALED UNITS BETWEEN PRINT OUTS 
I PRINT =0 OFF, I PR I NT = 1 ON 
HMIN=1 . 0E-9 
HMAX=H/100. 

HUSE=HM I N* 1 000 . 

I ERR=0 

INITIAL CONDITIONS 
Z 0=0.0 
YO ( 1 ) = 1 . 0 
Z 1 =0 . 0 
DG 9 1=2,6 
YG ( I ) =0. 0 

INITIALIZE FLOW, DENSITY, PRESSURE AND TEMPERATURE AT Z=0 
CALL FLOW ( ZO, TEMP, PRESS, FLOWR, DENSITY, DWDZ) 

3UESS AT INITIAL CONDITIONS FOR X(8> 

X70=P I N*W / (EPSNU*C) 

X70 IS RHO+ (0) WHICH IS GIVEN 

INITIAL CONDITION FOR RHO-(O) IS UNKNOWN — ASSUME COO VALUE 
X K 0 = X 7 0 * C 0 0 
Xe=XK0/X70 
YO (7) =X70/B 

IF ( I PRINT .EQ. 0) GO TO 300 
CALL FUN (N, ZO, YO, F) 

TEMP=TZZ+273. 

IF <Z0. EQ. 0. 0) CALL PRINT ( ZO, TEMP) 

WRITE (6, 191) 

FORMAT (///,T7, 1 HZ , T20, 4HCRI 1 , T32, 4HCR3 , T45, 5HCR2I , T57 , 4HC 121 , 

1 T69, 4HC I*I , T80, 4HC 1 1 , T91 , 6HCRHO+3 , T103, 6HCRH0-I , T1 12, 

2 9H INVERSION ) 

CONTINUE 

DO 10 1=1,7 
X ( I ) =B*YO ( I ) 

XS=XK0/X (7) 

X9=X <5>-,5*X (6) 

X7STAR=X (7) +X8 

USE SUBROUTINE SIGMA TO CALCULATE CROSS SECTION SIGMA 
CALL SIGMA (SIG2) 

IF ( I PR I NT .EQ. 0) GO TO 222 
WRITE (6, 199) ZO, (X (I) , 1=1,7) , XB, X9 
WRITE <6, 303) ZZZ,TZZ,PZZ, ETAZ Z , WZZ 
CALL FUN (N, ZO, YO, F) 

TEMP=TZZ+273. 

IF ( ZO. EQ. 1.0) CALL PRINT ( ZO, TEMP) 

CONTINUE 

SAVE THE DATA FOR FUTURE PLOT ROUTINES 

DATA ARRAYS ARE DATA (50, 50) AND FDPTDAT (50, 40) 
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IF ( I PRINT . EQ. 0) GO TO 227 
ICOL= (NG- 1 ) *6 
ND=ND+ 1 

DATA (ND, IC0L+2) =X ( 2) 

DATA (ND, ICOL+ 1 ) =Z0 
DATA (ND, ICOL+3) = X <4> 

DATA (ND, ICOL+4) =X (5) 

DATA (ND, ICOL+5) = X (6) 

DATA (ND, ICOL+6) = X9 
NNN= (NG-1 ) *5 
FDPTDAT (ND, NNN+1 > =Z0 


NNN+2) =WZZ 
, NNN+3 ) =ETAZ Z 


FDPTDAT (ND 
WRITE Z , R , 


FDPTDAT (ND 
FDPTDAT (ND 
FDPTDAT (ND, NNN+4) =PZZ 


X9 


6773 


303 


3567 

199 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

100 


ETAZZ, WZZ 
' X , T 1 5 


;HT= , F7. 3, 2X, T30, 


NNN+5 ) =TZ Z 
!2, I*, I, INV 

WRITE (8, 6773) ZO, X (2) , X (4) , X (5) , X (6) , 
WRITE (8,6773) ZO, WZZ , ETAZZ , F'ZZ , TZZ 
FORMAT ( 1 X , F6 . 3, 2X , 5 ( 2X , E 1 5 . 6 j ) 

CONTINUE 

IF (ZO.LE.O.O) GO TO 3567 
IF ( I PR I NT .EQ. 0) GO TO 223 
WR I TE ( 6 , 303 >ZZZ,TZZ,PZZ 
CONTINUE 

FORMAT <1X,T2,3HZ= ,F10.3, 

L 7HF'T0RR= ,F9.4,2X, 

l T55, 9HDENSITY = , F9. 6, 2X , TSO, 3HW= ,E14.7 
CONTINUE 

FORMAT ( 1 X , E12. 5, 8E12. 5, El 2. 5 ,E12.5 ) 


USE 7TH ORDER RUNGE KUTTA INTEGRATION SCHEME WITH VARIABLE STEP 
STEP SIZE CAN VARY FROM HMIN TO HMAX 
ZO IS STARTING VALUE FOR Z 

Z1 IS NEXT STOPPING POINT IN INTEGRATION SCHEME 
TOL IS TOLERENCE 

I ERR IS ERROR CODE TO DETERMINE IF INTEGRATION WAS SUCCESSFUL 


CONTINUE 
Z 1 =Z 1 +H 

IF ( Z 1 .GT. 1.00) GO TO 111 

CALL RKF7 (N, ZO, Z 1 , YO, TOL , FUN , PD , MTH , HM I N , HMAX , HUSE , WK , I ERR ) 
IF < I ERR .NE. 0) WRITE (6, 444) IERR, ZO, Z 1 , (YO ( I ) , 1 = 1 , 7) 

IF ( IERR .NE. 0) STOP 1717 


444 

FORMAT (IX, 18HIERR 

IS N 

200 

CONTINUE 
X<7)=B*Y0<7) 
X8=XK0/X (7) 


60 

FORMAT (1X,5(2X,E14 
GO TO 300 

.6) ) 

500 

CONTINUE 
DO 110 1=1,7 


110 

X ( I ) =B*YQ < I ) 
X8=XK0/X (7) 



,15, (4E16.7,/) 


X9=X (5) 5*X (6) 
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111 CONTINUE 

XX7L=X (7) 

C X (7) IS CALCULATED RHO+(L) VALUE 

XKCAL=R2*X (7) *X (7) 

XKCAL IS CALCULATED VALUE FOR XKO 

ABC IS DIFFERENCE BETWEEN CALCULATED AND INITIAL VALUE 
DIF= ( (XKO-XKCAL) /XKCAL) *100. 

ABC=DIF 

IF ( I PRINT .EQ. 0) GO TO 224 
WR I TE ( 6 , 202 ) D I F , X KCAL , X KO , COO 
24 CONTINUE 

02 FORMAT < 1 X, 13HD I FFERENCE = , El 8 . 9 , 2X , 1 2HXKCAL = , E1S.9, 

1 2X, 1 OH XKO = , E1S. 9, 2X , 6HC00 = ,E18.9 ) 

37 CONTINUE 

RHOM I N=RHO- (0) =XKO/RHO+ <0) =XK0/X70 
RHOM I N= X KO / X 7 0 
F OUT =RHOM I N*W*EPSNU*C 
GAINDB=10. *AL0G10 (POUT /PIN) 

GAIN=PQUT /PIN 

IF ( I PRINT .EQ. 0) GO TO 226 
WRITE (6, 193) R1 , R2, F', PIN, POUT, GAINDB, GAIN 

226 CONTINUE 

193 FORMAT ( IX, 5HR1 = , F 1 0 . 7 , 2X , 5HR2 = ,F10.7,1X, 

1 8HPRESS = , F 1 0 . 7 , 2 X , 6HF' I N = , E 1 4 . 7 , 3 X , 7HP0UT = ,E14.7,/, 

2 1 X , 9HGAINDB = , E 1 4 . 7, 3X , 7HGA IN = ,E14.7 ) 

501 CONTINUE 

NDMAX=ND 

C NDMAX IS THE MAXIMUM NUMBER OF DATA POINTS 

RETURN 
END 
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APPENDIX B 


OUTPUT FROM SIMULATION MODEL 



INPUT PARAMETERS 


PTO - 15.0 
OMEG1 - 3000 
CON - 1689 
COO - 2 , 13398E8 
R1 - 0.00 
R2 - 0.9584792 
LC - 1000 
ZOL - 500 
XNRHO - 1.0 
TO - 300.0 

OUTPUT 


P in 

P _ 
out 

Gain 

Gain db 

10.0 

128.4570 

12.8457 

11.0875 

5.0 

124.9000 

24.9801 

13.97595 

2.0 

122.2330 

61.1166 

17.8616 

1.0 

119.9611 

119.9611 

20.7904 

0.5 

116.9325 

233.8649 

23.6895 

0.2 

111.3596 

556.798 

27.4570 

0.1 

105.7950 

1057.954 

30.2446 

0.01 

80.3968 

8039.678 

39.0524 
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SCRLE DISTANCE, 


Z ildOJNI 



SCALE DISTANCE, 




DISTANCE, 




AXIAL DISTANCE, SCALED 
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